FUNCTION CA_IIR() provides an estimate of a filter's coefficients using a
stochastic gradient method.
Based on the paper "A complex adaptive algorithm for IIR filtering",
IEEE Transactions on Acoustics, Speech and Signal Processing, vol 34, no 5, 1986.
INPUT:
x: filter input [(N+1) x 1]
d: desired response [(N+1) x 1]
M: number of taps
mu: step-size
OUTPUT:
a, b : estimated taps [M x 1]
e: estimation error for each tap [1 x n]
y: filter output
y(n)= dot(a,y) + dot(b,x);
Complex Valued Nonlinear Adaptive Filtering toolbox for MATLAB
Supplementary to the book:
"Complex Valued Nonlinear Adaptive Filters: Noncircularity, Widely Linear and Neural Models"
by Danilo P. Mandic and Vanessa Su Lee Goh
(c) Copyright Danilo P. Mandic 2009
http://www.commsp.ee.ic.ac.uk/~mandic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You can obtain a copy of the GNU General Public License from
http://www.gnu.org/copyleft/gpl.html or by writing to
Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
...........................................0001 % FUNCTION CA_IIR() provides an estimate of a filter's coefficients using a 0002 % stochastic gradient method. 0003 % 0004 % Based on the paper "A complex adaptive algorithm for IIR filtering", 0005 %IEEE Transactions on Acoustics, Speech and Signal Processing, vol 34, no 5, 1986. 0006 % 0007 % INPUT: 0008 % x: filter input [(N+1) x 1] 0009 % d: desired response [(N+1) x 1] 0010 % M: number of taps 0011 % mu: step-size 0012 % 0013 % OUTPUT: 0014 % a, b : estimated taps [M x 1] 0015 % e: estimation error for each tap [1 x n] 0016 % y: filter output 0017 % y(n)= dot(a,y) + dot(b,x); 0018 % 0019 % 0020 % Complex Valued Nonlinear Adaptive Filtering toolbox for MATLAB 0021 % Supplementary to the book: 0022 % 0023 % "Complex Valued Nonlinear Adaptive Filters: Noncircularity, Widely Linear and Neural Models" 0024 % by Danilo P. Mandic and Vanessa Su Lee Goh 0025 % 0026 % (c) Copyright Danilo P. Mandic 2009 0027 % http://www.commsp.ee.ic.ac.uk/~mandic 0028 % 0029 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0030 % This program is free software; you can redistribute it and/or modify 0031 % it under the terms of the GNU General Public License as published by 0032 % the Free Software Foundation; either version 2 of the License, or 0033 % (at your option) any later version. 0034 % 0035 % This program is distributed in the hope that it will be useful, 0036 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0037 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0038 % GNU General Public License for more details. 0039 % 0040 % You can obtain a copy of the GNU General Public License from 0041 % http://www.gnu.org/copyleft/gpl.html or by writing to 0042 % Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA. 0043 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0044 % ........................................... 0045 function [a,b,e,y] = CA_IIR(x,d,M,mu) 0046 0047 L=length(x); % length of data 0048 N=L-1 ; 0049 StepA=1; % one-step ahead prediction 0050 0051 for stepAhead=1:StepA 0052 0053 y = (zeros(1,L) + i*zeros(1,L))'; 0054 e=y; 0055 0056 a = transpose(zeros(1,M-1) + i*zeros(1,M-1)); 0057 b = transpose(zeros(1,M) + i*zeros(1,M)); 0058 0059 psia= transpose(zeros(1,M-1) + i*zeros(1,M-1)); 0060 psib = transpose(zeros(1,M) + i*zeros(1,M)); 0061 0062 phia= transpose(zeros(1,M-1) + i*zeros(1,M-1)); 0063 phib = transpose(zeros(1,M) + i*zeros(1,M)); 0064 0065 0066 %begin of algorithm 0067 for n = M : N 0068 if (n+1+stepAhead)>N 0069 break 0070 end 0071 0072 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0073 % INPUTS 0074 u = x(n:-1:n-M+1) ; 0075 conju=conj(u); 0076 0077 if n==M 0078 y(1:M) = u(1:M); 0079 end 0080 0081 yy = y(n-1:-1:n-M+1) ; 0082 conjyy = conj(yy); 0083 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0084 0085 0086 y(n+stepAhead)= dot(a,yy)+ dot(b,u) ; 0087 0088 0089 e(n) = d(n+stepAhead) - y(n+stepAhead) ; 0090 conje = conj(e(n)); 0091 0092 %%%%%%%%%%%%Updates of sensitivities or gradients%%%%%%%%%%%%%%%%%%%%%%%%%% 0093 if n==M 0094 phia = conjyy; 0095 phib = conju; 0096 0097 else 0098 phia = [conjyy(1) + dot(conj(a),phia) ; phia(1:end-1)] ; 0099 phib = [conju(1) + dot(conj(a),phib(1:end-1)) ; phib(1:end-1)]; 0100 0101 0102 end 0103 0104 0105 %%%%%%%%%%%%Updates of coeffs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0106 if n< 101 % only 100 first samples to train, then use same coefficients 0107 mua = mu; 0108 mub = 5*mu; 0109 else 0110 mua=0; 0111 mub=0; 0112 end 0113 0114 a = a + mua*( e(n)*phia ); 0115 b = b + mub*( e(n)*phib ); 0116 end 0117 0118 end 0119 0120 0121 %Plotting of graphs 0122 figure,subplot(211), plot(real(y),'k-.') 0123 hold on,plot(real(x),'r') 0124 legend('CA IIR','Actual'), grid 0125 axis([1 length(x) min(real(x)) max(real(x))]), xlabel('Time (samples)') 0126 subplot(212), plot(imag(y),'k-.') 0127 hold on,plot(imag(x),'r'), grid 0128 axis([1 length(x) min(imag(x)) max(imag(x))]), xlabel('Time (samples)') 0129 0130